

METHOD AND APPARATUS FOR MODELING THREE-DIMENSIONAL 
ELECTROMAGNETIC SCATTERING FROM ARBITRARILY SHAPED THREE- 
DIMENSIONAL OBJECTS 


This invention was made with U.S. Government support from the Defense Advanced 
Research Projects Agency, under Contract No. MDA972-95-C-0021. The U.S. 
Government has certain rights in this invention. 


This application claims the benefit of priority to provisional application number 
60/213,102, filed in the United States on June 21, 2000, and titled "Well-conditioned 
Boundary Integral Equations and Discretization for Three-Dimensional Electromagnetic 
Scattering". 


The present invention relates generally to computer aided design and object modeling. 
More specifically, it relates to a computer program product, apparatus, and method for 
determining equivalent surface sources (or actual surface currents in the case of a metal 
20 object) on the surface of a closed, arbitrarily shaped object (e.g., a component of a 
vehicle such as an airplane part) that result from the imposition of an arbitrary time- 
harmonic incident field (e.g., a radar wave). 
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BACKGROUND 


Software utilizing boundary integral equation techniques to solve electromagnetic 
scattering problems (e.g., using the method of moments) has existed at least since the 
1960s. Such techniques have the potential for making very accurate radar cross-section 
predictions. However, practical application of the techniques to realistic problems is 
severely limited by poor scalability of computational resource requirements (i.e., CPU 
time and memory) relative to the problem size and complexity. With the relatively recent 
introduction of computational techniques involving high-order and "fast" mathematical 
methods (such as the fast multipole method (FMM)), the situation has changed 
drastically. Now, instead of the resource requirements for a problem described by N 
unknowns scaling in 0(N 3 ) time, they scale in 0(NlogN) time, much like the Fast Fourier 
Transform (FFT) algorithm. Consequently, it has become possible to handle a much 
wider array of problems of practical interest. 

Although the methods used have improved drastically, fast methods must still be used in 
conjunction with iterative solvers. If the set of linear equations that arises from 
discretizing a boundary integral equation is ill conditioned, the number of iterations 
required to obtain an accurate solution is not well controlled, and may be very large. The 
overall efficiency of the solution process suffers accordingly. In addition, a high degree 
of numerical precision is generally required in the solution process, causing the need for 
rounding, and resulting in a loss of solution accuracy. 
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Unfortunately, the most widely used integral equation formulation for closed, perfect 
electrical conductors (PEC), namely, the combined field integral (CFIE) equation, is not, 
in general, well conditioned. The ill conditioning stems from the fact that the CFIE 
necessarily includes an integral operator component that is "hypersingular". This 
5 component is, in fact, not simply an integral operator at all, but rather the gradient of an 
integral operator. Whereas, ordinary integral operators tend to filter out the numerical 
"noise" associated with a discretization, differential operators tend to amplify it. 

The eigenvalue spectrum of a hypersingular integral operator is unbounded, which leads 

O 

10 directly to an ill-conditioned linear system. The ill effects become evident in the form of 

m 

jy an increase in the required number of iterations as one refines the discretization of a 

j J target. Fine discretizations are often required to spatially resolve source variations or 

si geometric detail on the subwavelength scale. This difficulty with the CFIE (as well as 

5 

in with some other boundary integral formulations in computational electromagnetics) is 

i£f 15 known as the "low frequency problem." 

U 
M 

Various methods have been employed to overcome this problem. The most commonly 
practiced method is to multiply the CFIE matrix by a "preconditioner" matrix at every 
step of the iterative solver procedure. An effective preconditioner improves the condition 
20 number of a linear system, thereby reducing the number of iterations required for 

solution. Examples of preconditioners that have been utilized are the block-diagonal 
preconditioner, the incomplete LU preconditioner, and the sparse approximate inverse 
preconditioner. Preconditioners such as these are general purpose - in other words, they 
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can be applied with varying degrees of success to any matrix. Unfortunately, however, 
none of them takes specific advantage of the analytical properties of the CFIE, so none of 
them is ideally suited for use with the CFIE. 


5 Instead of preconditioning the equation, it is also possible to "post-condition" it. This 
involves defining the unknown part of the original linear equation to be the driving term 
for some new linear system. If the combined linear system is properly conditioned, then 
the problem is effectively solved because, once the new auxiliary unknowns have been 
solved for, it is straightforward to determine the solution to the original linear system. 

k 0 10 The moment method approaches that use "loop-star" basis functions fall into this 

W 

jfg category. While some success has been achieved with this technique, it is cumbersome 

iB 

ill with limited extensibility, and has been applied only to low-order, moment method 

ii discretizations. 

O 

in 

W 1 5 Therefore, it is desirable to provide a preconditioner method that is specifically tailored 


for use with the CFIE. More specifically, the preconditioner should be designed to 
minimize the number of iterations required for a solution. This is especially desirable 
when the use of a fine discretization to resolve source variations or geometric detail on a 
subwavelength scale results in an ill-conditioned linear equation (the "low-frequency" 
20 problem). Most specifically, it is desirable to provide a computer program product, an 
apparatus, and a method for modeling the electromagnetic response of a closed arbitrarily 
shaped three-dimensional object to an arbitrary time-harmonic incident field by means of 
a well-conditioned boundary integral equation (BIE) by discretizing the well-conditioned 
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BIE to provide and solve a well-conditioned finite dimensional linear system to 
determine the distribution of equivalent surface sources on the object. 
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SUMMARY OF THE INVENTION 
The present invention provides an improved method for modeling electromagnetic 
scattering from an arbitrarily shaped three-dimensional object to an arbitrary time- 
harmonic incident field. It allows for the determination of the scattering with an 
]0 10 improved degree of numerical accuracy, while also providing for increased 
computational efficiency. 


OB 

m 
m 


Vi ~ In a first embodiment, the present invention provides a computer program product, used 
Tl in a computer system having a processor, for modeling an electromagnetic response of an 

t 

i«0 15 arbitrarily shaped three-dimensional object to an arbitrary time-harmonic incident field 

Q 

by means of a well-conditioned boundary integral equation (BIE). The computer 
program product comprises recording medium; mans, recorded on the recording medium 
for producing a discretized representation of a well-conditioned BIE to provide a well- 
conditioned, finite-dimensional linear system whose solution approximates a distribution 
20 of equivalent surface sources on a model caused by the arbitrary time-harmonic incident 
field; and means, recorded on the recording medium, for solving the linear system to 
determine values of the equivalent surface sources; whereby solution of the well- 
conditioned finite-dimensional linear system provides improved computational efficiency 
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and solution accuracy. Additionally, this embodiment of the present invention may 
comprise means, recorded on the recording medium for modeling electromagnetic 
scattering from the three-dimensional object due to the time-harmonic incident field as 
determined from the equivalent surface sources. Preferably, the computer program 
5 product is tailored for determining the electromagnetic scattering from a time-harmonic 
incident field in the form of a radar wavefront and where the three-dimensional object 
represents a component of a vehicle. 

In another, more specific embodiment, the BIE produced by the computer program 
product is in the form of: 

(E excitation term)+/(H excitation term) = {(E integral operator term)+^(H integral operation term)}j 
where, y is a complex constant and J represents the unknown surface source distribution. 

In yet another, still more specific embodiment, the BIE produced by the computer 
program product is in the form of: 

kT s (ik)n x E inc + iZT a (ik)n • H inc = {/ [- T\ik) ° T T {k) + T\ik) o T L (k) - k 2 T s (ik) o T s (k)]- K + (k)}j 

In a further embodiment, the present invention provides an apparatus for modeling the 

20 electromagnetic response of an arbitrarily shaped three-dimensional object to an arbitrary 

time-harmonic incident field by means of a well-conditioned boundary integral equation 

(BIE). The apparatus comprises a computer system having a processor, a memory 

connected with the processor, and an input/output device connected with the processor 

for receiving parameters representing an arbitrarily shaped three-dimensional object and 
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an arbitrary time-harmonic incident field; means, operating in the processor and the 
memory of the computer for utilizing the parameters representing the arbitrarily shaped 
three-dimensional object and the arbitrary time-harmonic incident field for producing a 
discretized representation of a well-conditioned BIE to provide a well-conditioned finite 
5 dimensional linear system whose solution approximates a distribution of equivalent 
surface sources on a model caused by the time-harmonic incident field; and means, 
operating in the processor and the memory of the computer for solving the linear system 
to determine values of the equivalent surface sources; whereby solution of the well- 
conditioned finite dimensional linear system provides improved computational efficiency 

O 

*fl 10 and solution accuracy. Similar to the computer program product embodiments, the 
00 

P apparatus further comprises means, operating in the processor and the memory of the 

jU 

J 3 ? computer for modeling electromagnetic scattering from the three-dimensional object due 

to the time-harmonic incident field as determined from the equivalent surface sources. 

jj? Again, the apparatus is preferably tailored for determining the electromagnetic scattering 
!-* 

Ifl 1 5 from a time-harmonic incident field represented in the form of a radar wavefront and 
U 

I "* where the three-dimensional object represents a component of a vehicle. 

In another, more specific embodiment, the BIE produced by the apparatus is in the form 
of: 

20 (E excitation term)+ /(H excitation term) = {(E integral operator term)+ y(H integral operation term)}j 
where, y is a complex constant and J represents the unknown surface source distribution. 
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In yet another, still more specific embodiment, the BIE produced by the apparatus 
product is in the form of: 

kT s (ik)n x E inc + iZT\ik)n ■ H " c = (f [- r B (*) o r r (Jfc) + 7*(iJfc) o r'ot) - Jfc 2 r* (*) o 

The present invention further provides a method for modeling an electromagnetic 
response of an arbitrarily shaped three-dimensional object to an arbitrary time-harmonic 
incident field by means of a well-conditioned boundary integral equation (BIE). The 
method operates on a computer system having a processor, a memory connected with the 
processor, and an input/output device connected with the processor for receiving 
parameters representing an arbitrarily shaped three-dimensional object and an arbitrary 
time-harmonic incident field. The method comprises the steps of: utilizing the 
parameters representing the arbitrarily shaped three-dimensional object and the arbitrary 
time-harmonic incident field for producing a discretized representation of a well- 
conditioned BIE to provide a well-conditioned finite dimensional linear system whose 
solution approximates a distribution of equivalent surface sources on a model caused by 
the arbitrary time-harmonic incident field; and solving the linear system to determine 
values of the equivalent surface sources; whereby solution of the well-conditioned finite- 
dimensional linear system provides improved computational efficiency and solution 
accuracy. 

Additionally, the method can further comprise the step of modeling electromagnetic 
scattering from the three-dimensional object due to the time-harmonic incident field as 
determined from the equivalent surface sources. The method is preferably operated 
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where the time-harmonic incident field is a radar wavefront and wherein the three- 
dimensional object represents a component of a vehicle, in order to determine the 
scattering of radar waves from a vehicle component due to in incident radar wave. 

5 In another, more specific embodiment, the BIE produced by the method is in the form of: 

(E excitation term)+ /(H excitation term) = {(E integral operator term)+ ^(H integral operation term)}j 

where, y is a complex constant and J represents the unknown surface source distribution. 

-^10 In yet another, still more specific embodiment, the BIE produced by the method is in the 
(39 

jjg form of: 

m ^(/*)nxE' nc +/zr^ 

5 DESCRIPTION OF THE FIGURES 

JT-5 

1 5 FIG. 1 is a block diagram representing a general purpose computer system used in 

□ 

u conjunction with the present invention; 

FIG. 2(a) is an illustrative diagram showing an object having a diameter of 2.5A, (where X 
represents a wavelength), and corresponding to element 300 of FIG. 3; 
FIG. 2(b) is an illustrative diagram showing the object of FIG. 2(a) modeled as a plurality 
20 of triangular patches, and corresponding to element 302 of FIG. 3; 

FIG. 2(c) is an illustrative diagram showing the object of FIG. 2(b) with an incident field, 
and corresponding to element 304 of FIG. 3; 
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FIG. 2(d) is an illustrative diagram showing the object of FIG. 2(c) with an incident field 
and with electromagnetic fields resulting from equivalent surface sources due to the 
incident field, and corresponding to elements 306, 3 08, and 310 of FIG. 3; 
FIG. 3 is a flow diagram demonstrating the method of the present invention in 
conjunction with the steps of inputting object and incident field parameters, dividing the 
model of the object into a plurality of patches, and the optional step of determining 
electromagnetic fields; and 

FIG. 4 is an illustrative diagram of a computer program product according to the present 
invention, stored on CD or DVD. 


DETAILED DESCRIPTION 
The present invention provides a computer program product, a method, and an apparatus 
for generating and solving well-conditioned boundary integral equations for modeling 
jfi three dimensional electromagnetic scattering of an incident wavefront from an object. It 

\H 1 5 should be understood that this description is not intended to limit the invention. On the 

0 

contrary, the invention is intended to cover alternatives, modifications and equivalents, 
which are included in the spirit and scope of the invention as defined by the appended 
claims. Furthermore, in the following detailed description of the present invention 
numerous specific details are set forth in order to provide a thorough understanding of the 
20 present invention. However, it will be obvious to one of ordinary skill in the art that the 
present invention may be practiced without the specific details. 


Introduction 
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In the presence of an incident wavefront, the surface properties of an object affect 
scattering of the wave from the object's surface. By properly modeling the object and the 
incoming wavefront, it is possible to determine how waves will be scattered from the 
object's surface. Determination of this scattering is of particular importance with regard 
to designing objects having reduced radar cross-sections as well as in designing antennae 
or other objects that are configured to maximize electromagnetic interception. 

The present invention, in a general sense, provides a computer program product, a 
method, and an apparatus for receiving a computer model of an object comprised of a 
mesh of a plurality of patches; for specifying an incident wavefront having an excitation 
wavelength, wavefront geometry (e.g., plane wave, etc.), and incident angle; and, based 
on the properties of the object and the wavefront, determining the equivalent surface 
sources resulting from the interaction of the wavefront with the object, and from the 
equivalent surface sources, modeling the scattering of energy from the surface of the 
object. 

Invention Overview 

The following overview is presented with reference to the figures in order to provide a 
high-level conceptual understanding of the present invention prior to discussing the 
mathematical details and specific advantages inherent therein. First, an overview of the 
data processing system of the present invention is provided, followed by a brief 
illustration and flow diagram explaining the operation of the present invention. 
Subsequent to the overview, further details of the present invention will be provided, 
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followed by a set of numerical results demonstrating the advantages of the modified 
CFIE of the present invention over the conventional CFIE. Finally, an appendix and 
glossary are provided to give the reader a convenient source for understanding some of 
the terms and conventions utilized herein. 

Data Processing System 
The present invention is designed to be operated on a data processing system, which may 
comprise a general or special purpose computer. Having reference now to the drawings, 
in FIG. 1, there is shown a computer or data processing system generally designated by 
reference character 100. The computer 100 generally includes a central processing unit 
(CPU) 102 which is coupled with a memory 104. The memory 104 generally comprises 
both random access memory (RAM) and read-only memory (ROM) portions. The 
computer 100 further comprises at least one input/output device 106 connected with the 
CPU 102, such as a keyboard, a pointing device, a scanning device, or other available and 
useful devices. A display 108 is also connected with the CPU 102 to provide output to a 
user. Further, a storage mechanism 110 is provided to enable a user to store and retrieve 
data. Non-limiting examples of various storage mechanisms 110 that may be of use with 
the present invention include hard drives, optical drives (such as a compact disc or digital 
versatile disc), and magnetic drives (such as floppy disks, magnetic tapes, and Zip™ 
disks). Finally, a communications adapter 112 (such as a modem or network interface 
card) may also coupled with the computer 100 in order to allow communication with 
other computers 100. Note further that the data processing system may also comprise a 
network of computers either local to or remote from each other. In addition, the various 
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components of the data processing system may all be local or may be distributed 
remotely from each other. 

Although described at a major-component level, it is assumed that the computer 100 
5 includes all of the necessary adapters to drive all of the peripherals necessary for use with 
the present invention. Additionally, various commercially available processors could be 
used in conjunction with the present invention, and could be programmed to operate the 
logic of the preferred embodiments. 

Illustration and Flow Diagram of the Operation of the Present Invention 
Referring to FIG. 2, and FIG. 3, the operation of the present invention is described by 
means of an illustrative description and a flowchart. After starting operation of the 
present invention, represented by the begin 300 element in FIG. 3, an object 200 is 
divided into a plurality of patches, represented by element 202 in FIG. 2 and element 302 
in FIG. 3. Note that methods for dividing the object into a plurality of patches are well- 
known in the literature, and that the patches may be of an arbitrary geometric shape (e.g., 
triangles, squares, etc.). Next, parameters are supplied for the incident field and the 
properties of the object as depicted by element 304 of FIG. 3. The incident field is 
depicted as element 204 in FIG. 2. The dashed lines surrounding elements 302 and 304 
of FIG. 3 are intended to indicate that they are environmental to the present invention, as 
the patches 202 representing the object 200 and the parameters for the incident field 204 
and the object 200 may simply be passed as inputs to a system incorporating the present 
invention. Optionally, however, the present invention may include a system for dividing 
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the object 200 into a plurality of patches 202 and for more directly determining the 
associated parameters. In addition, although the object 200 is represented as a sphere and 
the incident field 204 is represented as a series of planar waves, these specifics are 
provided for illustrative convenience only. The object 200 may be any closed arbitrary 
5 shape and the incident field 204 may be any arbitrary time-harmonic field. Next, as 
shown by element 306 in FIG. 3, the system develops a discretized representation of a 
well-conditioned boundary integral equation (which takes the form of a well-conditioned 
matrix equation), modeling the interaction between the incident field 204 and the object 
200. The well-conditioned linear system is then solved numerically to determine 

o 

\Q 10 equivalent surface sources on the object, as depicted by element 308 in FIG. 3. After 

5 

\M solving for the equivalent sources, the present invention may utilize the equivalent 

W sources in order to determine electromagnetic fields resulting therefrom, as depicted by 

^ element 310 of FIG. 3. The dash-dot-dot line in element 310 is intended to indicate that 

I! 

P 

1^ it is optional. The electromagnetic scattering from the object is represented by element 
IH 15 206 in FIG. 2. 

(CSS, 
U 

\ + 

By providing a well-conditioned linear system in the determination of the equivalent 
sources, the present invention helps to ensure that the solution to the linear system can be 
performed with fewer iterations and with increased numerical precision relative to 
20 conventional methods. In fact, without the use of the well-conditioned linear system, it is 
sometimes impractical to determine the equivalent surface sources accurately due to the 
limited numerical precision of one's computer. 
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Recording Medium 

In one embodiment, the present invention takes the form of a computer program product 
recorded on a recording medium. A disc, shown in FIG. 4 represents a recording 
medium 400 such as a compact disc (CD) or a digital versatile disc (DVD). Although 
5 depicted as a CD or DVD, any recording medium could be utilized, examples of which 
include magnetic discs and tapes, and other optical storage medium. 


Further Details of the Present Invention 

Well-Conditioned CFIE Formulation 

Q 

\Q 10 In accordance with features of the present invention, to provide improved computational 

IS 

W efficiency and accuracy, modifications are made to the combined field integral equation 
jU 

'™ (CFIE) for the solution of electromagnetic scattering problems in three dimensions. 

*«■! 

J* Unlike the conventional CFIE, the modified CFIE of the present invention is well- 
jfj conditioned. It combines, a standard magnetic field integral operator with a pre- 

jjj 15 conditioned electric field integral operator, with the resulting modified CFIE being an 

□ 

M integral equation of the second-kind that has no spurious resonances. As a result of the 
pre-conditioning, the number of iterations required to solve the CFIE stabilizes even as 
discretization of the surface is refined, whereas refinement of the discretization in 
conjunction with the conventional CFIE results in an ill-conditioned linear system and a 
20 growing number of iterations. The conventional CFIE, not being well-conditioned, 
suffers from two major drawbacks. 
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As mentioned previously, there are three principal reasons why well-conditioned 
boundary integral formulations are beneficial. The first is the fact that sparse 
decompositions of the integral operators of scattering theory, but not their inverses, exist, 
so iterative solvers must be employed. Well-conditioned systems of equations can be 
solved with few iterations. Second, using a fine discretization to resolve source 
variations or geometric detail on a sub-wavelength scale typically results in an ill- 
conditioned linear equation as well (also known as the "low-frequency" problem in 
computational electromagnetics). Third, because the process of generating a solution to 
the set of discretized equations can require a very wide range of values (in terms of 
magnitude), the resulting solutions may be quite inaccurate because a high degree of 
numerical precision may be required, which may exceed the accuracy permitted by the 
data processing system. 

Only second-kind integral equations, or objects with similar spectral behavior (such as 
appropriately preconditioned differential equations) can be solved with fully controlled 
approximation error. The correct operators are the sum of a constant (or at least well- 
conditioned and easily invertible) operator and a compact operator. 

Boundary integral operators of scattering typically violate this requirement in one of three 
ways: the spectrum may accumulate at zero; the operator may have an unbounded 
spectrum; or the operator may have small eigenvalues associated with resonances, often 
unphysical (e.g., "spurious resonances"). 
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For electromagnetic scattering from perfectly conducting (PEC) surfaces, the standard 
boundary integral equations are the electric field integral equation (EFIE) 

- n x E inc = TJ (where inc stands for incident), (1) 

and the magnetic field integral equation (MFIE) 


ZnxH mc = 


IB 

IB 10 where the integral operators Tand K are defined by 

fu 

m n= T(k)j 


= ikn(x)x J^'jG(*;x,x'^ and (3) 


\ / I I \ 7 7 / \ / j 

3 s I * 


s 

m 


jH KJ = K(k)J = - n(x) x J<ft' VG(k; x, x') x J(x') , (4) 

£i S 
W 

|== 

where V denotes differentiation with respect to x, and n(x) is the unit normal to the 
1 5 surface at x. Note that the other terms follow the usual conventions: J = Zn x H is the 
unknown surface current, E mc and H mc are the incident electric and magnetic fields, 


respectively, Z = ^Jifs is the wave impedance, and G{k\ x, x') = exp(/fcr)/ 4/zr is the 3d 
Helmholtz kernel with r = |x - x'| being the distance separating field and source points. 
Harmonic time dependence e ~ 16 * is assumed. 


20 
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The MFIE is a second-kind integral equation. Unfortunately, this equation is suitable for 
an unacceptably small class of electromagnetic problems. It is inapplicable to open 
surfaces, becomes ill-conditioned in the presence of geometric singularities, and suffers 
from spurious resonances. The EFIE has both a compact piece and a hypersingular piece 
(coming from the double gradient term). The spurious resonances of the MFIE may be 
eliminated by adding the EFIE to form a combined field integral equation (CFIE). The 
cost of doing so is the introduction of EFIE's hypersingular piece, which spoils the 
conditioning for fine discretizations (or low frequencies). 

It is known that a hypersingular integral operator and a first-kind integral operator are 
ideal preconditioners for each other, in the sense that the composition of the two has the 
spectral characteristics of a second-kind integral operator. A similar approach may be 
employed to analytically precondition the EIFE. In fact, as disclosed herein, the electric 
field integral operator T preconditions itself. 

Two issues are important for the successful application of this idea to closed bodies. 
First, only the local (or short distance) behavior of the preconditioner is important for 
asymptotic conditioning. Thus, the EFIE may be preconditioned by multiplying it by an 
electric field integral operator corresponding to an arbitrary wavenumber, real or 
complex; if the wavenumber has a positive imaginary part, the introduction of any 
additional resonances may be avoided. If the EFIE preconditioner reproduced the MFIE 
resonances, then the CFIE would also have them. Second, one must take care that the 
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discretization of the product of the preconditioner and preconditioned operators preserves 
the correct spectral properties. 

The present invention provides well-conditioned formulations for closed surfaces. 
5 Example numerical results are provided further below in order to demonstrate the 
advantages of the modified CFIE formulation over the conventional CFIE. 

In the literature, integral operators have been constructed from the kernel for the Laplace 
and Helmholtz equations in 2d, Through the construction of these operators, it has been 

O 

\Q 1 0 observed that the product of a first-kind operator, constructed from an undifferentiated 

m 

\M kernel, and a hypersingular operator, constructed from a twice differentiated kernel, has 

Hi 


the desirable spectral characteristics of a second-kind operator. Since the EFIE integral 

|«* 

s 

operator T has both of these, it may be seen that the composition of two such operators 

m 

its? 

jjj T 2 = T o T would include a constant operator and a compact operator. Also of concern is 

jjl 1 5 that the product of hypersingular components would produce another hypersingular 
P 

i a operator. However, it is easy to see that the rotation operation nx in the definition of T 9 

as shown in equation (3) above, annihilates the component of the surface vector field 
normal to the surface and ensures that the product of the two hypersingular operators is 
identically equal to zero. Indeed, applying the hypersingular component of the second T 
20 operator to an arbitrary tangential surface vector function f(x') produces a surface 
gradient function, 

n x V0(x) = -(nx V) \ds* VG(Jfc; x, x') • f (x') , (5) 
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(note the dot-product operator, •) which the hypersingular component of the first T 
operator, in turn, annihilates (for closed surfaces) by virtue of the identity 

V s -[nxV#x)] = 0, (6) 
5 with V s denoting the surface gradient operator on S; the identity shown in equation (6), 
the surface analog of the 3d identity V • [V x <j>{x)\ = 0 , can be found and is valid for any 
sufficiently smooth function ^on S. It follows immediately from equations (3), (5), and 
(6) that T 2 behaves as a second-kind integral operator. 


O 


1^ 10 In this specification, the spectral properties of the EFIE and MFIE integral operators and 
W 

^ combinations thereof for the PEC sphere, a simple 3d target for which the spectral 
153 
\± 

'si 

p discussion may be expanded to incorporate other object shapes. The use of the PEC 


properties of these operators are known analytically, are discussed. However, this 


M sphere herein is intended to provide a non-limiting example of an object to which the 

m 

O 15 present invention may be applied. A complete set of basis functions on the surface of a 
sphere of radius a is given by the vector spherical harmonics, 

Xfc.(*,0- . ; * "xVyjM and (7) 

/V/(/+i) 

V lm (O,<p) = nxX lm (0,<p), (8) 
defined here in terms of the scalar spherical harmonics Y lm (0,(p) . 

20 The result of applying T and K + =(K + ^) to each basis function is: 
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ra)IM-!" J ' (fc,)H ' (fa0U 4 (9) 

Kj k (ka)«, (fa>)X,J 


and 


fxj f Q\ (ka) M, (ka) X /m 


™\v:nimimij- (10) 

where J, and H / are Riccati-Bessel and (first-kind) Riccati-Hankel functions of order /, 
and k is the wavenumber associated with the kernel of each integral operator. The 
Riccati-Bessel and Riccati-Hankel functions are defined in terms of spherical Bessel and 
Hankel functions j)(x) and h\ x \x) by 

J/W^OOand (11) 
U { {x) = xff\x) (12) 


Although chosen basis functions X lm and U/ OT are not eigenfunctions of the operator T{k\ 
they are eigenfunctions of T 2 (k) = T(k) o T(k) : 

T2{k tfv'"} = - J ' {ka) m < {ka) j; {ka) m; (Ab) {u| m } • (13) 

The operator T*(k) has a bounded spectrum, since, in the limit of large /, its eigenvalues 
accumulate at -1/4 (a result which follows from the asymptotic properties of ;*/ and h\ X) ). 
However, as is evident, from equations (10) and (13), the operator T 2 {k) also shares 
resonances (at the zeros of J\ (ka) for the X /m modes, and at the zeros of J, (ka) and for 

the U /m modes) with the MFIE operator K + (k) . This is also evident from the identity 
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T *( k ) = I _ K 2 (k) = -K~(k) o K + (k), (14) 

where K~ = — - K . Therefore, although T 2 (£) is a second-kind integral operator, it is 
2 

not a suitable component of a resonance-free combined field integral equation for closed 
bodies. 


As previously mentioned, boundedness of the spectrum of the product of two EFIE 
operators (of the form shown by equation (3)) is assured if they have the same short- 
distance behavior, a condition that does not require the two operators to share the same 
\Q wavenumber (propagation constant). If EFIE operators are chosen with different 

m 

^ 10 wavenumbers, T(t) and T(k 2 ) , it is possible to simultaneously obtain a bounded 


ii 


product and avoid MFIE resonances. 


IS The following analysis indicates that ik is a particularly good choice for the wavenumber 
in the preconditioning operator (assuming that the wavenumber k is real). The 


□ 

l== 1 5 eigensystem for T (ik) o T(k) on a sphere is 


r m (r AM \ f > {ika) H ' {ika) 3 ' {ka) m,{ka) x * V (15^ 


It is straightforward to show (given the properties of ji and hj l) ) that the eigenvalues of 
T (ik)oT(k) accumulate at ^ and - l - for the X /m and U /m eigenmodes, respectively, 
and that T (ik) o T(k) does not share any resonances with the MFIE operator K + (k) . 


20 


Page 24 of 55 

6/15/2001 
HRL053FINAL 


# • 


Since T (ik) o T(k) is a second-kind integral operator (in the sense described below) and 
does not share any resonances with K + (k), it is possible now to develop the well- 
conditioned CFIE operator. The simplest form of such an operator is 


5 T{ik)oT(k) + aK+{k), (16) 

where a is a constant to be chosen. In developing this modified CFIE operator, the EFIE 
portion has been modified before being added to the MFIE portion (which is already a 

second-kind integral operator). The same applies to the excitation side of the equation. 

P 

|| 10 The resulting CFIE may be written as 


m 

09 


w 

in 


- T(ik) (n x E' nc ) + aZn x H " c = [T(ik) ° T(k) + aK+(k)]j . (17) 


The eigensystem for the CFIE operator in (16) is 


% ^iuJ = -i [ j 1 (* 0 )H,( ifa )j;M +fa J,( fa ,))H;( fel )uJ- (18) 

15 

If one chooses a = ±1 then, as a function of the argument ka, these eigenvalues have no 
zeros. For a = +1 , they circle the origin of the complex plane. 

Other well-conditioned (modified) CFIE operators can be devised, for example, by 
20 preconditioning the MFIE part before combining it with the preconditioned EFIE part. 
Two example forms of alternative well-conditioned (modified) CFIE operators are as 
follows: 
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T(ik) o T(k) + aK + (ik) o K + (k) 
and 

T(ik) o T(k) + anxK + (ik) onx K + (k) . 


(19) 


(20) 


It has been observed that the numerical behavior of the three modified CFIE operators (in 
(16), (19), and (20)) is similar. 

By virtue of the analysis on spheres, the modified CFIE operators in (16), (19), and (20) 
have been shown to be second-kind and resonance-free. However, given that the 
asymptotic behavior of the eigenvalues on a smooth surface stems from the short distance 
behavior of the kernel, it may be seen that the asymptotic behavior of the various 
operators on spheres should also obtain for any closed surface that can be obtained by 
smooth deformation of a sphere. The numerical results presented further below support 
this argument. Results are also presented for a cube, which, like many targets of practical 
interest, has geometric singularities. These suggest that the modified CFIE formulations 
should be well-conditioned for a wide class of closed surfaces. 

Discretization and an Alternate Form of the Modified CFIE Formulation 

Discretization of the modified (well-conditioned) CFIE formulation utilized in the present 

invention may take place by any appropriate method. Suitable methods are common in 

the literature, non-limiting examples of which include the Nystrom discretization method 

and the method of moments. A particularly useful reference in regard to the Nystrom 

discretization method, which may be used in conjunction with the present invention may 
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be found in reference [13], in the Background, titled "Numerical solution of the 
Helmholtz equation in 2d and 3d using a high-order Nystrom discretization" by Canino et 
al. Note also, incorporating the approach of the present invention into existing "fast" 
solvers is completely straightforward. 

There are several ways to discretize the product operator T (k x ) o T(k 2 ) . The simplest 
and most straightforward approach, multiplying the discretized representations of the 
individual operators, can lead to numerical difficulties. The reason is that it is relatively 
difficult to make the discretized representations of the hypersingular part of each operator 
sufficiently accurate (especially for high-spatial-frequency eigenmodes) to numerically 
effect the cancellation that obtains analytically. 

Effective discretizations of T (k x ) o T(k 2 ) can be obtained either by discretizing the 
product operator directly or by reformulating the product operator to eliminate the 
product of hypersingular operators. The second approach has been implemented using a 
reformulated product operator that eliminates all instances of hypersingular operators. A 
short derivation of the reformulated equation is provided below. 

The first step toward obtaining a more useful form of the product operator T o T{k 2 ) 
is to separate each integral operator into its singular and hypersingular components. 
Introducing the abbreviations, 


T x =T(k x ) and 


(21) 


T 2 =T(k 2 ), 


(22) 
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and substituting, the following can be written: 



1 ' 


(23) 



2 » 


(24) 


where, 


T*J = n(x) x jds'G(k m ; x, x')J(x'), and 


(25) 


S 


CJ S (n(x) x V) \ds' VG(*„; x, x') ■ J(x') . 


(26) 


5 


The product operator 7] o T 2 can be expanded into four terms. Each of the two cross 
terms, T } s o T 2 H and T X H o7 2 s , can be transformed (by Stokes's theorem) into a product of 
new, single-gradient integral operators on S plus a line integral around the boundary of 5. 
The term formed by the product of the hypersingular integral operators, T X H o T 2 H , reduces 
to a line integral. The result is further simplified by noticing that two of the three line 
integrals, when applied to J, can be combined into a single term whose argument is 
identical to the incident electric field E mc by virtue of equation (1). 

The next step is to reformulate the excitation side of the equation, taking advantage of the 
fact that the incident wave obeys Maxwell's equations. By applying Stokes's theorem, 
the term T } H [n(x') x E l>IC (x ! )] may be rewritten as the sum of a single-gradient integral 
operator on V x E wc (x') and a line integral that exactly cancels the line integral 
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m 


involving E inc on the other side of the equation. A further simplification follows from 
Faraday's Law, V x E = ico/jK . 

The final result for the analytically preconditioned EFIE with reformulated integral 
5 operator product is 

- iktfin x E') - Z^TT(n ■ H') = o T 2 T + ^Tf o T 2 L - k x k 2 T x s o 7* - ^ o J 

k 2 k x k 2 fc 2 

(27) 

where the various integral operators are defined by: 

Ty = n(x)x \ds<S/G(k m , x, xW)> (28) 


IS 

iU 10 Tfy = n(x) x [ds'n(x')x V'G(k m , x, x'^(x') , (29) 

IS s 


fi 

o 
□ 


T m L f = VG(k m ; x, x') ■ f(x'), (30) 

S 

T T J = n(x) • \ds'VG(k m , x, x 1 ) x f(x') , (31) 

s 

r m s f ^ n(x) x j"<fc'G(fc m ; x, x'XOO , and (32) 

r.V s n(x) x Jd'G(* m ; x, x')tfx') , and (33) 

as 

1 5 where m=l , 2. Note that T° , T% , and T* map scalar functions to surface vector 
functions, whereas T£ andfj do the reverse. The operator on the right hand side of 
equation (27) maps surface vector functions into surface vector functions. 


Page 29 of 55 

6/15/2001 
HRL053FINAL 


Now, closed surfaces will be discussed, and it may be observed that T x oT 2 behaves like 
a second-kind integral operator. 

If S is a closed surface, the term T X E o T 2 J vanishes and equation (27) simplifies to 

- iktfQi x E') - Z^-T x a (n • H') = S(k x , k 2 )J , (34) 
k 2 

where 

5, 2 ^ S(k x , * 2 ) = o Tl + ^Tf o - °r. (35) 
Now, several features of Sn will be discussed. 

First, all of the individual operators that comprise S 12 involve kernels with one or no 
gradients on the Helmholtz Green's function G. All such integral operators are bounded. 

Second, the eigenvalues of the integral operator Sn do not accumulate at the origin. This 
may be demonstrated by examining its three components, T" o T 2 , Tf o T 2 , and 

o T 2 . The operator T x a o T 2 is a second-kind integral operator for the transverse 
(divergence-free) component of J, and is identically zero for the longitudinal 
(irrotational) component of J. Likewise, the operator Tf ° T 2 is a second-kind operator 
for the longitudinal component of J, and is identically zero for the transverse component 
of J. Since any surface current distribution can be decomposed into longitudinal and 

k k 

transverse components, the sum — T x a oT 2 +—Tf oT 2 is a second-kind integral 

k x k 2 
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operator; subtracting k x k 2 T x s ol^a compact operator, does not change this result. As 
observed previously, resonance sharing may be avoided by setting k x = ik and k 2 =k. In 

this case, the eigenvalues of S12 accumulate at two points, ± — , rather than only at - . 


5 Third, the spectrum of Si 2, after discretization, is bounded and includes accumulation 
points at the expected locations. However, an accurate discretization will have zero (or 
very small) eigenvalues wherever the EFIE operator T(k 2 ) has a resonance. Thus, it has 
to be combined with an appropriate discretization of the MFIE operator, to obtain an 

0 

k Q effective discretization of the CFIE. 

in 

m 10 

m 

W Finally, it should be noted that equation (34) is manifestly insusceptible to the "low- 

^ frequency" problem that plagues the EFIE. Since the well-conditioned behavior of Sn 

it 

□ k k 

0 1 comes from the composite operators —T x o T 2 and -^Tf o T 2 L , both of whose 

I** k x k 2 

m 
o 


prefactors have modulus unity (assuming \k^ = |& 2 | = k ), and since the term 

1 5 k x k 2 T { s o T 2 tends to zero as k -> 0 , the full operator Sn remains well conditioned in the 
limit of low frequency. 

In summary, although the operators T x o T 2 and Sn have identical spectral properties for 

closed bodies, it is easier to construct an accurate discretization for S12 because it is 

20 composed of less singular integral operators. Matrix representations of S/2 have bounded 

spectra, but also suffer from spurious resonances inherited from the EFIE operator T(k 2 ) . 
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These resonances can be eliminated by combining Sj 2 with K + (k 2 ) (or the modified 
MFIE operators in equations (19) and (20)). The result is a well-conditioned system of 
linear algebraic equations. 

An Alternate General Formulation of the Modified CFIE 
The boundary integral equation used for the present invention may be formulated more 
generally, as follows: 

(E excitation term)+/(H excitation term)= {(E integral operator term)+/(H integral operation term)}j 
where, 

y is a complex constant, 

J represents the unknown surface source distribution, and 

the (E excitation term) and (E integral operator term) are of a formulation selected from 
the group consisting of A and B, with 

A representing the combination of 

(E excitation term) = T{k x )(- n(x) x E inc (x)) and 

(E integral operator term) = T(k x ) T(k 2 ) , and with 
B representing the combination of 


(E excitation term) = 



(E integral operator term) = 


Page 32 of 55 

6/15/2001 
HRL053 FINAL 


^T\k,)oT T (k 2 ) + ^T\k 2 )oT L (k 2 )-^ 

the (H excitation term) and (H integral operation term) are of a formulation 
selected from the group consisting of A', B', and C\ with 
A' representing the combination of 

(H excitation term) = Zn(x) x H mc (x) and 

(H integral operation term) = K + (k 2 ); 
B' representing the combination of 

(H excitation term) = K\k x )(Zn(x) x H' nc (x)) and 

(H integral operation term) = K + (k x )o K + (k 2 ); and 
C representing the combination of 

(H excitation term)= n x K + (k x )(Zn(x) x H' wc (x)) and 

(H integral operation term) = [n x ^(fc^Jo [n x £ + (£ 2 )], wherein; 
T(k) = ikn(x) x jds^G(k; x, x')J(x) + -p- v[VG(£; x, x') • J(x')]| , 

K + (k) = ~ + K(k),and 

K(k) = -n(x) x fa VG(£; x 3 x') x J(x') , and 

s 

wherein J = Zn x H is the unknown surface current; 
£ is a complex number; 

E" c and H mc are the incident electric and magnetic fields, respectively; 

Z = yjju/z ls ^ wave impedance; 
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G{k\ x, x 1 ) = exp(ikr)/4nr is the 3d Helmholtz kernel, 

where r = |x - x'| is the distance separating field and source points with 

harmonic time dependence e~ 16 * assumed, and where 

T a (k)# = n(x)x jdsWG(k; x, x'#(x') , 

s 

T p {k)(f> = n(x)x J<fc'n(x')x V*G(i; x, x'^(x') , 

5 

r L (*)f s Jds'VG(*; x 5 x f )f(x'), 

s 

T T (k)f == n(x) • VG(*; x, x') x f (x') , and 

s 

T s (k)f = n(x) x jds'G(k; x, x')f (x f ) , and where 

s 

V and V denote differention with respect to the field point x and source 
point x' , respectively, and n(x)and n(x')are outward-pointing unit 
vectors normal to the target surface S at x and x' , respectively, and 

k 2 is the propagation constant in the medium external to S; 
k { is a complex number chosen to avoid resonances. 


This general form may be modified to provide the preferred form of the modified CFIE, 
as previously mentioned, 

kT s (ik)n x E inc + iZT a (ik)n - H ,nc = jf [- T a (ik) o T T (k) + 7*0*) o T L (k) - k 2 T s {ik) o T s (k)]- K + (k)}j 
, where (ik) is substituted for (k) in the various parts of the equation as defined above. 

Specifically, for example, the term 
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T a (k)</> = n(x) x jds } VG(£; x, x'^(x') is simply changed to 

s 

T a (ik)</> = n(x)x jds'VG(ik; x, x'^(x') . 

s 

The same type of change is applied to all of the other terms. 

5 Numerical Results Demonstrating the Advantages of the Modified CFIE Over the 
Conventional CFIE 

Below, the computational performance of the conventional CFIE (referred to below as 
CFIE1), 

- n x (n x E') + Zn x H' =[nx T(k) + K + (£)Jj (36) 


is} 10 is compared with the preconditioned CFIE (referred to below as CFIE2) 

U kT s (ik)n x E' + iZT a (ik)n ■ H' - Zn x H' 

'-4 

" = {/ [-T a (ik)oT T (k) + T fi (ik)oT L (k)-k 2 T s (ik)oT s (k)\-K + (k)}j (37) 


produced by combining equation (17) (with a = -1) and equation (34) (with k x = ik and 
k 2 = k ). Individual operators in these equations were discretized using a high-order 
1 5 Nystrom scheme. In all cases, the wave impedance Z was set to unity. 

Three examples are provided below to demonstrate the operation of the present invention. 
The first example shows how the condition number of each operator, defined as the ratio 
of the largest to smallest singular values, depends on the fineness of discretization. Table 
20 1, below, lists the condition number (CN) of the matrix representing each CFIE operator 
as the size of the sphere decreases. In all cases, the same discretization was used, created 
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by placing a 6-point quadrature rule on each of the 80 nearly identical patches that cover 
the sphere, for a total of 960 unknowns. As the sphere radius decreases, the condition 
number for the CFIE2 integral operator stabilizes at about 2, whereas the condition 
number for the CFIE1 continues to grow in inverse proportion to the radius. 


radius(X) 

CFIE1 

CFIE2 

1 

4.2 

3.04 

1/4 

15 

2.68 

1/16 

59 

2.04 

1/64 

230 

1.99 

1/256 

940 

1.97 

1/1024 

3800 

1.97 

1/4096 

15000 

1.97 


5 

»0 Table 1 : Condition number of CFIE matrices for shrinking PEC spheres 

!» 

IB 

9 

|2 The second test compares iterative solver performance for the modified CFIE with that 
for the conventional CFIE. The target geometry consists of two PEC spheres, one with a 

P 

iJI 10 radius of XI 2 , the other set at a resonant radius, namely, the first zero of JJ (2nrlX) or 

|- 

iH r ~ 0.43667457A . The spheres are separated by a X l\ 00 gap. The patches near the gap 

P 

,ss were subdivided by a factor of 10 to adequately resolve the currents, which vary rapidly 
there. Table 2 compares the iteration counts and radar cross section (RCS) errors for 
several discretizations. The iterations columns list the maximum number of iterations a 
1 5 conjugate gradient squared (CGS) route required to reach a residual error of 1 0" 3 . A 

solution computed from a substantially more refined discretization provided an accuracy 
reference. The stated error is the root mean squared (RMS) value of the difference 
between the monostatic RCS of the comparison solution and the reference solution at 

181 angles. For identical discretizations, the two methods had about the same error. The 
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data show a dramatic difference, however, in the iteration count behavior of the two 
methods in response to discretization refinements. 


unknowns 

patches 

CPIE1 

CPIE2 



iterations 

error 

iterations 

error 

1496 

748 

60 

0.46 

9 

0.40 

4488 

748 

126 

0.18 

11 

0.18 

996 

498 

44 

0.61 

9 

0.48 

2988 

498 

103 

0.23 

12 

0.16 

5976 

498 

163 

0.016 

11 

0.029 


Table 2: Iteration count and solution error vs. discretization for two PEC spheres. 

The third test also compared iterative solver performance for the two CFIE formulations. 
In this case, the target was a cube of size U (where X is the wavelength of the radiation). 
Numerical results are presented for five different discretizations, the first of which was 
derived from a mesh (i.e., a set of patches) obtained by dividing each face into four 
squares. The second mesh was constructed from the first by subdividing each square into 
four smaller squares. The third mesh was constructed from the second by subdividing 
edge-touching patches in half along a line parallel to the edge; patches adjacent to two 
edges (i.e., corner patches) were divided into quarters. Meshes for the fourth and fifth 
discretizations were constructed by recursively applying the procedure by which the third 
mesh was constructed from the second. This process, known as patch tapering, is useful 
for resolving the source singularities that arise in the vicinity of geometric singularities. 
It also puts stress on the conventional CFIE because points near the edges get close 
together. Table 3 lists the maximum and average number of iterations the CGS routine 

required to obtain solutions for 92 independent excitations to a residual error of 10" 3 . 
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The total number of unknowns is the result of using a 9-point quadrature rule on each 
square or rectangular patch. The iteration count for CFIE2 grows very slowly with 
increasing taper depth, whereas for CFIE1 it increases steadily as expected. 


unknowns 

taper 

CPIE1 

CPIE2 


depth 

max 

ave 

max 

ave 

432 

0 

12 

6.5 

10 

4.3 

1728 

1 

18 

9.9 

11 

4.9 

3888 

2 

26 

14 

11 

4.9 

6912 

3 

41 

23 

11 

5.5 

10800 

4 

58 

36 

13 

5.7 


IS 
139 

m 


i ' 

in 


15 


20 


Table 3: Iteration count vs. taper depth for U PEC cube. 
Appendix - Second-Kind Integral Operators 

The standard definition of a second-kind integral operator is an operator of the form 
M + K, (38) 
where X is a constant, / is the identity, and K is a compact operator. In scattering theory, 
one encounters operators of the form 

Vi+A^+tf, (39) 
where \ and ^ are constants and P\ and P 2 are orthogonal projection operators such that 
P X +P 2 =I. (40) 
Operators of the form shown in equation (39) possess most of the desirable properties of 
second-kind integral operators. Although somewhat non-standard, such expressions are 
referred to as second-kind operators throughout this description. 
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Glossary of Terms and Conventions 

Composition Operator - The composition operator used herein, and represented by the 
symbol, o , describes the following operation. Given a first integral operator, 

A = jds'K A (x 9 x') , and a second integral operator, B = jds"K B (x', x") , then the 

s s 

composed operator A o B = jds'K A (x, x') jds"K B (x', x") . 

s s 

Equivalent surface sources - According to a well-known mathematical theorem (the 
"equivalence principle"), the fields outside a volume V can be completely specified by a 
distribution of equivalent sources on the boundary of the volume S. If S corresponds to 
the surface of a metal target (which approximates a perfect electrical conductor), the 
equivalent sources correspond to true, physical surface currents. 

Means - Generally, for purposes of this invention, the term refers to a computer program 
or sub-program, or a portion of a computer program operative for performing the task 
that the means is for. The means is generally considered to be language and computer 
platform independent, and in the form of computer instructions (or code) that may be 
either in the form of software or hardware. 

Object - This term object refers to primarily any arbitrarily shaped closed surface for 
which equivalent surface sources and the resultant electromagnetic scattering are to be 
determined. Objects of interest could include, as non-limiting examples, aircraft or parts 
thereof, automobiles, tanks, boats, missiles, and fixed structures. The equations set forth 
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herein are developed for perfectly electrically conducting (PEC) materials or 
approximately perfectly conducting materials such as metals. However, extension to 
surface scattering with more general boundary conditions is considered routine and 
within the scope of the present invention. The extension to an interface between two or 
more dielectrics, for example, is straightforward, with the resulting operators having 
condition numbers slightly lower than the case described herein. Although primarily 
developed for closed surfaces, the approach can potentially be extended for the case of 
open, perfectly conducting surfaces. 
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